#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#define TIME 500
#define OUTPATH "ssaout.txt"
#define INDMAX 70
#define STEPCOUNT 180000	
#define SPECIES 46
#define REACTIONS 70


/*
	Switching mechanism;
	swState :
		0 ... no switch;
		1 ... MCT(iEC) -> BFP(iPC)
		2 ... BFP(iPC) -> MCT(iEC)
*/
int swState = 1;
int swTime	= 90000;


unsigned z1, z2, z3, z4;  


/*
unsigned TausStep(unsigned z, int S1, int S2, int S3, unsigned M)  
{  
  unsigned b=(    (unsigned)pow((z << S1), z)     >> S2);  
  return z = pow(((z & M) << S3) , b);  
} 
*/ 

unsigned TausStep(unsigned *z, int S1, int S2, int S3, unsigned M)  
{  
  unsigned b=(((*z << S1) ^ *z) >> S2);  
  return *z =  (((*z & M) << S3) ^ b);
}  

unsigned LCGStep(unsigned *z, unsigned A, unsigned C)  
{  
	return (A*(*z)+C);  
}  

float HybridTaus()  
{  
  return 2.3283064365387e-10 * (              // Periods  
    TausStep(&z1, 13, 19, 12, 4294967294UL) ^  // p1=2^31-1  
    TausStep(&z2, 2, 25, 4, 4294967288UL) ^    // p2=2^30-1  
    TausStep(&z3, 3, 11, 17, 4294967280UL) ^   // p3=2^28-1  
    LCGStep(&z4, 1664525, 1013904223UL)        // p4=2^32  
   );  
} 

int main(int argc, char *argv[]){

	FILE * fp;
	fp = fopen(argv[1],"w");

  //Switching
  swState = atoi(argv[3]);

	/*Parameters*/
	
	
	//Load vectors
	//float w[REACTIONS]		= {0};
	//float c[REACTIONS]		= {0};
	//int species[SPECIES];
	//#include "ssaspecies.h"
	//#include "ssacoefs.h"
	#include "model.h"

	//Seed random
	srand(time(0));

	//High log
	#define HL_AGM 20
	#define HL_PI 3.14159
	#define HL_LN2 0.69314
	#define HL_MBITS 8
	float hl_f;
	float hl_a;
	float hl_g;
	float hl_at;
	float hl_s;
	float hl_x;
	float hl_y;
	float hl_lnf;

	//Misc
	int i		= 0;
	int j		= 0;
	int mi		= 0;
	int tstep	= 0;
	float tau	= 0;
	float t	= 0;
	
	//Ssa misc
	float r1;
	float r2;
	float wsum = 0;
	float wthresh[REACTIONS];
	float wtcumsum = 0;

	
	//State change vector; (Read only memory)
	//#include "ssadif.h"
	

	//HybridTaus
	z1 = 131; z2 = 129; z3 = 130; z4 = rand();
	
	/*
	FILE * fpx;
	fpx = fopen("ssadif.h", "w");
	for (i = 0; i<REACTIONS; i++){
		fprintf(fpx, "{");
		for (j = 0; j < SPECIES; j++){

		  	dif[i][j] = -reactants[i][j]+products[i][j];
			if(j!=SPECIES-1) fprintf(fpx," %d,", dif[i][j]);
			else fprintf(fpx," %d", dif[i][j]);
		}
		if(i!=SPECIES-1) fprintf(fpx, "},");
		else fprintf(fpx, "}");
		fprintf(fpx,"\n");
	}

	fclose	(fpx);
	t = 0;
	return;
	*/
	
	
	for (i = 0; i<REACTIONS; i++){for (j = 0; j<SPECIES; j++) dif[i][j] = -reactants[i][j]+products[i][j];}
	
	
	
	while(tstep < STEPCOUNT){

	//r1 = ((float)rand())/RAND_MAX;
		//r2 = ((float)rand())/RAND_MAX;
		r1 = HybridTaus();
		r2 = HybridTaus();
	
		//Include propensities calculations
		#include "ssaweights.h"
	
		//Calculate wsum
		//weight thresholds
		wsum=0;
		wtcumsum = 0;
		for (i = 0; i<REACTIONS; i++) 	wsum = wsum + w[i];
		if (wsum <= 0){break;}
		for (i = 0; i<REACTIONS; i++){
			wtcumsum = wtcumsum + w[i]/wsum;
			wthresh[i] = wtcumsum;
		}


		//Determine tau;
		//tau = 1.0/wsum * logf(1.0/(1-r1));

		//High log
		hl_f = 1/(1-r1);
		hl_s = hl_f*pow(2,HL_MBITS);		
		hl_y = 4.0/hl_s;
		hl_x = 1;
		//agm start
		hl_a = 0.5*(hl_x+hl_y);
		hl_g = sqrt(hl_x*hl_y);
		for(j=0; j<HL_AGM; j++){
			hl_at = hl_a;
			hl_a = 0.5*(hl_a+hl_g);
			hl_g = sqrt(hl_g*hl_at);
		}
		//agm end
		hl_lnf = HL_PI / (2*hl_a) - HL_MBITS * HL_LN2;

		//Determine tau.
		tau = 1/wsum * hl_lnf;

		

		//Determine reaction index.
		mi = -1;
		for (i = 0; i<REACTIONS; i++)	mi = (r2 < wthresh[i] && mi == -1)*i + (r2 >= wthresh[i] || mi != -1)*mi;

		//Update species matrix;
		for (i=0; i<SPECIES; i++) species[i] = species[i] + dif[mi][i];


		//Print to file
		//fprintf(fp,"%f\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%f\t%d\t%f\t%d\t%d\t%f\t%f\t%f\t%f\n",t, species[BFP], species[mCITRINE],species[P3],species[P4],species[TALA_KRAB],species[TALB_VP16],species[iPC],w[20], mi, w[mi],species[TALA_VP16], species[TALB_VP16], r1, r2, w[3],w[13]);
		//fprintf(fp,"%f\t      %d\t%d\t%d\t           %d\t%d\t%d\t       %d\n",t, species[BFP], species[TALA_VP16], species[TALB_KRAB],   species[mCITRINE], species[TALB_VP16], species[TALA_KRAB],  mi);
		for (i=0; i<SPECIES; i++) fprintf(fp,"%d ", species[i]);
		fprintf(fp," %d\n", mi);

		t = t+tau;
		tstep++;
	
		//Switch mechanism;
		if(tstep == (int)STEPCOUNT/2){
			//1 ... MCT(iEC) -> BFP(iPC)
			if(swState == 1){
        //printf("switch 1! %d\n", tstep);
				species[iPC] = INDMAX;
				species[iEC] = 0;
				//species[E_KRAB] += species[iEC_E_KRAB];
				//species[iEC_E_KRAB]  = 0;

			}
			//2 ... BFP(iPC) -> MCT(iEC)
			if(swState == 2){
        //printf("switch 2! %d \n", tstep);								
				species[iPC] = 0;
				species[iEC] = INDMAX;
				//species[PIP_KRAB] += species[iPC_PIP_KRAB];
				//species[iPC_PIP_KRAB]  = 0;
			}
		}



	}

	
	//printf("%f, %f\n",t,tau);
	printf("%d, %d\n", species[BFP], species[mCITRINE]);
	fclose(fp);
	return 0;
}


